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ABSTRACT 


Parameter estimation algorithms are developed in the frequency domain for systems modeled 
by input/output ordinary differential equations. The approach is based on Shinbrot’s method 
of moment functionals utilizing Fourier based modulating functions. Assuming white meas- 
urement noises for linear multivariable system models, an adaptive weighted least squares 
algorithm is developed which approximates a maximum likelihood estimate and cannot be 
biased by unknown initial or boundary conditions in the data owing to a special property 
attending Shinbrot-type modulating functions. Application is made to perturbation equation 
modeling of the longitudinal and lateral dynamics of a high performance aircraft using flight- 
test data. Comparative studies are included which demonstrate potential advantages of the 
algorithm relative to some well established techniques for parameter identification. Deter- 
ministic least squares extensions of the approach are made to the frequency transfer function 
identification problem for linear systems and to the parameter identification problem for a 
class of nonlinear time-varying differential system models. 
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1. Introduction 

Estimating aircraft parameters from flight or wind-tunnel data has been an important 
modeling activity for aerodynamicists over several decades. Various kinds of parametrized 
models exist to encompass the steady and unsteady flight conditions (Klein, 1989). Here the 
focus is placed on estimating the parameters of ordinary differential equation input/output 
models for which a number of methods exist. These include: approximations like the 5 
operator, block pulse and orthogonal function expansion methods; moment functional tech- 
niques; filtering techniques like the maximum likelihood/extended Kalman filter, and the 
conversion of an identified discrete-time model (obtained by a variety of methods) into a 
continuous-time model using, for example, the bilinear transformation. The books by 
Unbehauen and Rao (1987) and Johansson (1993) are good introductions to the identification 
of continuous-time models, and Sinha and Rao (1991) contain chapters by researchers report- 
ing recent results on some of these methods. 

The method of moment functionals, also known as the modulating function technique, is 
characterized by the use of integration-by-parts on an equation error model of the system to 
transfer the continuous -time derivatives on the input/output variables to derivatives on a set of 
smooth user-chosen ‘modulating’ functions. The result of this process is a set of algebraic 
equation errors which is characterized by functionals on the data that have to be computed 
before initiating the parameter estimation. If the modulating functions are of the Shinbrot 
type (Shinbrot, 1954,1957), then the chosen functions must satisfy certain end point conditions 
for the time interval over which data is given and, as a consequence, the resulting equation 
error vector will be devoid of unknown initial or boundary conditions on the data. Herein lies 
one advantage of Shinbrot’s method in that it effectively decouples the state and parameter 
estimation problems for data collected under transient conditions. Utilizing Fourier based 
modulating functions, another advantage is that the system identification problem can be 
equivalently posed entirely in the frequency domain, and the functionals can be calculated by 
efficient DFT/FFT techniques. It is Shinbrot’s method that is promulgated in this report. 

Following a discussion in Section 2 of the aircraft models used to motivate the formula- 
tion and the main assumptions concerning the data, Shinbrot’s method is developed in Section 
3 for single input single output (SISO) linear differential systems based upon a chosen set of 
Fourier modulating functions. Using a second order example, a comparison is made via simu- 
lation with a commercially available prediction error algorithm which illustrates the improved 
accuracy accrued for the modulating function technique over a range of signal to noise ratios. 
The formulation of Section 3 for SISO models is extended in Section 4 to the multivariable 
(MIMO) aircraft models introduced in Section 2, and results are presented for a high perfor- 
mance aircraft based upon actual flight data. Using the same data sets, a comparison is made 
with the modeling results obtained via a seasoned time domain based maximum likelihood 
identifier. This comparison indicates good corraboration when both yield useful models, but it 
also supports the efficacy of the modulating function technique for data sets in which the time 
domain based identifier failed. The final section discusses various extensions. Included here 
is showing how Shinbrot’s method can be developed for a frequency transfer function 
identification problem, and an extension of the theory to a class of nonlinear input/output 
models. 



2. Aircraft Perturbation Models 

A model for longitudinal motion can be represented by the pair of linear differential 
operator equations: 

Aj(p)a(r) -B x {p)u{t) (1) 

A l (p)q(t) = B 2 (p)u(t) (2) 

where {A^ip), B ]{p), B 2 {p)) are polynomials in the differential operator p=d/dt whose 
coefficients represent parameters to be estimated given sampled versions of the input signal 
i/(0 and the two output signals (a(0, <7(0) over some time interval [0,7]. Here (a(r), q{t)) 
denote respectively the angle of attack and the pitch rate. In one application u(t) will be the 
horizontal tail deflection, denoted by Sy, (f ), while in another u (?) will be the longitudinal stick 
deflection denoted by T|/, (t). The aircraft is operating in a closed loop as indicated in Fig. 1. 
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Fig. 1. Longitudinal dynamics block diagram 

A model for lateral motion, presumed to be independent of the longitudinal motion, can 
be represented by the three differential operator equations: 

A 2 {p )(3(r ) = B u (p )b a (r) + B l2 <p ) 6 , (r) (3) 

A 2 (p)p(0 =B 21 (p)b a (t) + B22(p) 5 r(0 ( 4 ) 

A 2 {p )r (0 = B 31 (p ) 6 a (0 + B 32 (p ) 6 r (t ). (5) 

Here the problem is to estimate parameters comprised of the coefficients in the polynomials 
(A 2 (p), Bjj(p)), / =1,2,3 and j— 1,2, given sampled versions of the two input signals 
(b a (r), 5 r ([)), which are respectively the aileron and rudder deflections, and the three output 
signals ((3(f), P (0. f(f)). which are respectively the sideslip angle, rolling velocity and yaw- 
ing velocity, over a time interval [0,7], 

Underlying each set of models (l)-(2) and (3)-(5) is a state vector differential equation of 
the generic form: 

i(t) = Fx(t) + Gu(t), y(t) = Hx (t ). (6) 

In terms of transfer functions, the following model relation holds: 

-77 -B(s) = H(sl-FT 1 G (7) 

A (s) 
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where A (s )=det(s7 -F ) is designated by A rfs) for the longitudinal model (l)-(2), A 2 (s) for 
the lateral dynamics model (3)-(5), and the B ( 5 ) matrix is comprised of the polynomials £, ( s ) 
or ( 5 ) as appropriate to each set. The reason for working with the above particular model 
forms, i.e., the various outputs within each set are operated on by the same scalar-valued 
A {p) operator, is a result of the ease with which the modulating function technique developed 
in Pearson et al. (1993b,c) for single input single output (SISO) systems can be extended to 
this form. In each set of models, the polynomial A {p ) is chosen monic with order fixed on 
the basis of physical considerations: Specifically for the perturbation equations relative to 
steady flight, the orders of (A ^s), A 2 (s)) are chosen to be (2, 4) respectively, and the polyno- 
mials comprising B (s ) are chosen of order one less than the corresponding A (s ) polynomial. 
This means that the model (l)-(2) possesses a total of 6 coefficient parameters, while the 
model (3)-(5) possesses 28 total coefficient parameters. In some applications, the lateral 
dynamics has included a fourth output variable - the Euler roll angle <J >(/ ) - in which case the 
lateral dynamics possesses 36 total coefficient parameters. In addition, there are unknown 
weighting parameters that arise in the estimation procedure as will be indicated below. 

2.1. Data Assumptions and Bandwidth Considerations 

The major assumptions concerning the input/output data [n (t ),y(t )], 0<.t<.T , are: 

(i) oversampling is employed in terms of the Nyquist sampling theorem, i.e., the sampled 
data [u (kT/N),y ( kT/N )] , k =0, 1 • • N, is collected in sufficient quantity and rate that the fol- 
lowing inequality prevails: 

1F b < < F s = j (7) 

where F B is the bandwidth of the system, 

(ii) the various input signals are observed with negligible measurement noise and with 
sufficient frequency content to avoid degeneracy in the least squares estimates, and 

(iii) noises corrupting the output signals are modeled as additive zero-mean independent 
white Gaussian noises with unknown variances. In addition, all stochastic calculus operations 
are presumed to be carried out in the mean square sense. This validates the various modula- 
tion processes on the differential equation models, provided all the pertinent signals are mean 
square differentiable of order equivalent to the model order. 

The first two assumptions are quite reasonable for the aircraft modeling considered here, 
viz., the data is collected at a 50 Hz sampling rate, the system bandwidths are expected to be 
on the order of a few Hertz, and the input signals are relatively free of measurement noise 
while possessing sufficient frequency content for the estimation problem. The third assump- 
tion facilitates specifying a weighted-least-squares cost function in the frequency domain, via 
Fourier modulating functions of the Shinbrot type, whose minimization approximates a max- 
imum likelihood estimate. This formulation for SISO systems is summarized in the following 
section. The basic calculation needed to set up the estimation problem is computing a finite 
number (specifically M+n where the integers (M,n) are explained below) of Fourier series 
coefficients for each data variable. If /(r), (Kt <T, is any such data variable and 
f k =f(kT/N), *=0,1 • • N, its sampled version, then the following staircase approximation to 
the Fourier integral is presumed to be carried out in obtaining the needed coefficients F (m)\ 
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(8) 


where / =n^T, coQ=2n/7 , m=0, 1, • • M +n and W N = e~ ,2nlN . 

In the above, coq is the ‘resolving’ frequency which is a user-chosen parameter if the data 
length T is selectable, n is the order of the polynomial A (p ) specified for the model; M is 
an integer that can be chosen on the basis of the number of unknown parameters or on the 
basis of bandwidth considerations if r=27t/co 0 is fixed. Thus, if Hq denotes the number of 
unknown parameters (n e =6 for the model (l)-(2) and n e =28 for the model (3)-(5)), M can be 
chosen by the guideline: 

M ~ 2n e ~ 4 n Q (9a) 


N-l 

*«0 


which will provide about 2 n Q to An 0 algebraic equations in the discrete-frequency domain 
upon which to base the least squares estimate. If the data interval [0,7"] is selectable, then the 
guideline: 

(M +n )co 0 ~ co B => T = 2n/co 0 ~ 2n(M+n)lw B (9b) 


where © B =2 nF B is the system bandwidth. Hence, these two guidelines, which were upheld 

for the aircraft modeling carried out here, will uniquely determine the pair (A/.coq). 2 Note that 
if the number of required harmonics, M+n , is artificially extended to equal the number of 
available samples, N, then the right hand side of (8) is the discrete Fourier transform of the 
data sequence f k on the sampled interval, thus facilitating a DFT/FFT algorithm. Calculating 
the DFT’s of the input/output data is a negligible computation in all the applications encoun- 
tered thus far. 


In terms of Hertz, the left hand side of (9b) means (M+n )F 0 ~ F B , F q=HT , which 
together with the sampling assumption (7) implies the inequality: 


M_+_n 

N 12 


<< 1 . 


(9c) 


This inequality implies that only the lowest few percent of the available DFT harmonics will 
be utilized in the least squares algorithm. In particular, fewer than 10% were utilized in the 
aircraft modeling applications. Hence, high frequency noise rejection is inherent in the 
Fourier based modulating function technique if the guideline (9b) is upheld. It can be noted 
that this inequality is also consistent with a basic property of discrete Fourier transforms relat- 
ing to the first and last halves of the DFT being complex conjugates of one another which, in 
turn, leads to the requirement that (M+n) < N 12 in order to preserve uniqueness for the 
Fourier coefficients. Furthermore, inequality (9c) is consistent with a fundamental tenet of 
identification in that the number N of available samples will ordinarily greatly exceed the 


2 It is advantageous in the case of high order systems to normalize the [0,7" ] interval to [0,2n] so 
that ©(pi, or some value dose to unity. The reason is that derivatives in the time domain model are 
replaced by powers of coq in the frequency domain model. Calculating powers like ©o', will lead to 
numerically better conditioned regression equations in the frequency domain when © 0 ~1. This will be 
clear from the formulation in Section 3. 
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number n e of unknown parameters which from (9a) translates into the inequality N > >M +n . 


3. Shinbrot’s Method Via Fourier Modulating Functions 

The basic idea of Shinbrot’s method will be illustrated with respect to the problem of 
estimating the parameters for a second order single degree of freedom mechanical system. 
Given the force and displacement signals, represented respectively as input/output data 
IM0j(0] on a time interval [0,1], find parameters (a 1 ,a 2 ,h 1 ) that best fit (in some sense) 
the model: 


y(t) + a l y(t) + a^y(t) = b 1 u(t) i 0 <t<T. 

Multiplying the equation error for this model by a smooth function <|>(0 and integrating over 
the given data interval: 



(0 + <*jy(0 + tfiKO - biu(t) \dt = e 


where e represents the accumulated equation error over [0,7], Using integration-by-parts on 
the derivative portions of the above integral, the ‘modulated’ equation error s is equivalent to 
T T T T 


e = 




)y{t)dt - a if<K0y(r>* + a 2 \§{t)y(t)dt - b l \^{t)u{t)dt + BC 


0 


0 


0 


where BC stands for ‘boundary conditions’; specifically, 

BC = {<K0y(0 - 4>(0 y( 0 + <*i4>(f)y(0}J. 

As introduced by Shinbrot within the context of this second-order case, <J)(r ) is a modulating 
function, herein to be called Shinbrot-type in order to distinguish from other types, if it 
satisfies the four end-point conditions: 

«j)(0) =4»(7) =(j>(0) = <j>(7) = 0. 


Then BC= 0, i.e., presuming [u (t),y(t),y(t)] are bounded on [0,7], and hence these conditions 
imply that the modulated equation error reduces to 

T T T T 

-b x ^(t)u ( t)dt 



= Y 0 “fliYi -a 2 Y2 — ^iY3 

which involves the linear functionals y h /=0, 1,2,3, defined on the input/output data, but not 
derivatives of the data. Repeating this process with a sequence of linearly independent 
smooth functions (t ), each satisfying the same four end-point conditions, will lead to a set 
of algebraic equation errors e(m ), m=l, • ■ M, characterized by linear functionals y^m) on 
the data, which can be used as a basis for estimating the coefficients (a h a 2 ,b : ). This esti- 
mate cannot be biased by unknown initial or boundary conditions in the data owing to the 
special end-point conditions of the Shinbrot modulating functions. Other moment functional 
techniques, such as the Poisson functionals (Saha et al. 1991), do not share this property. 
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Although Shinbrot’s method was developed in an aeronautical setting (Shinbrot, 1954), it 
has remained relatively obscure in relation to the maximum likelihood and extended Kalman 
filtering techniques that have evolved for modem aerospace applications (Klein, 1993). There 
may be several reasons for this but in the main they all stem from a lack of clear guidance as 
to a good choice of Shinbrot-type modulating functions. Consider in this regard the following 
order-n Fourier modulating function of the Shinbrot-type which is defined and represented 

equivalently by 3 


* 


m 71 


(0 = l. e ~ ima *(e~ i<s>of -l) n , 0 <Lt<> 7 = 2ti/cd 0 


= 




c k = (-1) 


n-k 


0 


fe] 


i m+n # 

1 kmm 


(10a) 

(10b) 

(10c) 


where i='PT, CDo=2n/7 is the resolving frequency, and m is any integer which shall be 
referred to as the ‘modulating frequency index’. The family of functions defined by (10) for 
any specified index set me {m lt • • m M } possesses several properties, four of which are as 
follows: 


Property 1. For k =0, 1, • • n- 1: p*<J» w/I (0 = 0 at t = 0 and t = 7 

where p denotes the differential operator, i.e., p=d/dt, p=djdt, etc. 


( 11 ) 


Property 2. For a sufficiently smooth function z(t) defined on [0,7] and a differential opera- 
tor P{p) of order n , the integration of <j) OT Jt {t)P {p )z (t ) over [0,7] satisfies: 

T m+n 

(12a) 


, n {t)P{p)2 (t )dt = Y. C k-m P O'* ®o) Z (* ) 
0 kmm 


= A n P (im (£>q)Z (m) (12b) 

where Z(m) is the m th harmonic Fourier series coefficient of the function z (t ) on [0,7], and 
A" is the n ,h order finite difference operator, i.e., with Q (m )=P (im to 0 )Z (m ), 

k ■ m+n 


A n Q(m)= £(-1)* \" \Q(n+m-k)= Z c k-mQ(k). 

km 0 L J kmm 


Property 3. For a zero mean continuous-time Gaussian white noise process v (t ) with covari- 
ance Ev (t j)v (r 2 )=a 2 5(f 2 -r i) and a differential operator P(p) of order n, the complex- valued 
stochastic sequence e(m), m =0,1 • • , defined by 


3 Equation (10b) follows from (10a) by use of the binomial expansion where j denotes the bi- 
nomial coefficient, and (10c) follows from (10b) by a change in the summation index. 
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e(m) = ^ mjl (t)P{p)v{t)dt = e R (m ) + /e/(An) 


is also Gaussian with the covariance relations among its real and imaginary parts given by: 


Ee R (m)e R (m+l ) = — 


m+n 


(P(0))%[m]4o[/] + \l n [l] £ co 0 )| 

k~m+l 


E ej ( m )£/ (m+l) = 


a 


m+n 


-c P (0)) 2 p 0 [m ]p<)[/ ] + [/ ] £ c*_ m Cfc—fn \p (ik co 0 )l 

J:-m +/ 


E e R (m )e 7 (m +1 ) = 0 


where the notation }!„[/] is defined by the string of unit pulses: |l n [/ ] = 1 


1 for (Kl <n 
0 otherwise 


Thus, Hot" 1 ] and (l 0 [/] are unit pulses at the origin which corresponds to the zero frequency. 

Property 4. For sufficiently smooth functions z(t) and w(t) defined on [0,7’] and a 
differential operator Pip) of order n, the integration of § mn (t)w(t)P (p)z(t) over [0,7'] 
satisfies 

T 

.n (t ) w (' )P (p )z it )dt =W(m )® A w P (im co 0 )Z (m ) (13) 

where ® stands for linear convolution in the discrete frequency domain of Fourier series 
coefficients, i.e., with Q {m )=A" P {im (o 0 )Z (m ), 

Wimy&Qim) = £ W(l)Q(m-l). 

l*z-oo 

Remarks on Verification, (i) Property 1 follows directly from the representation (10a) 
since ie~ l<i>cf —l) n is the n ,h power of a unit vector centered at -1 in the complex plane which 
starts at the origin and makes one complete revolution as t ranges over [0,7’]. It is this pro- 
perty which qualifies „ (f) as a Shinbrot-type modulating function of order n . 

(ii) The second property follows using the representation (10b) or (10c) while employing 
integration-by-parts on the left side of (12) and taking into account (11); during the process of 
transferring the derivatives on z(t) to derivatives on the modulating function <J) m n (t ), all 
boundary point evaluations are zeroed such that 

T T 

[<iw(OP(p)z it)dt = 

In view of the representations (10b, c), coupled with the fact that P (—p)e~^ =P (k)e~^ , the 
right side of (12) is obtained which verifies Property 2. This property facilitates trading time 
derivatives on data-related functions with finite differences in the frequency domain unencum- 
bered by unspecified initial/boundary conditions on the data. 

(iii) The third property hinges on the presumption that all signals of interest in the 


jz(t)P(-p)ty m n (t)dt. 
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identification problem are mean square differentiable of order n, the order of the underlying 
model. This property is proven in Pearson and Shen (1993c) drawing upon mean square 
results from the stochastic calculus. 

(iv) Property 4 follows using the Fourier series representation for w (t ) in the left hand side of 
(13), interchanging the summation and integration order, then applying Property 2 to each 
term in the sum. Property 4 is consistent with the Fourier analysis property involving the 
multiplication of two functions in the time domain and its corresponding relation to convolu- 
tion in the frequency domain. 


3.1. Parameter Estimation for SISO Models 


Consider the single input single output system configuration in Fig. 2 where the transfer 

function H(s) is represented by a ratio of polynomials ~7~r- 

^4 (5 ) 


U(t) 


H(s) 





Fig. 2. System with disturbance v(r) 

The system is modeled on [0,7] by the differential operator equation: 

P n y( r) + E a jP n ~ J y(0 = £bjP n ~ J u( t) + e(t) ( 14 ) 

7- 1 7-1 

where e(t) represents the equation error. For consistency between Fig. 2 and the model, e(t) 
is defined via the relation: 

e(t) = A(p)v(t) = Y,ajp a - j v(t), a 0 = 1. ( 15 ) 

i - 0 

The problem is to estimate the (a ; ,bj ) parameters, j— 1,2 ■ n, which are coefficients in the 
polynomials (A (s)Ji(s)), given the data [u (t),y(t)] on [0,7] and assuming the output y(f) is 
corrupted by a zero mean white Gaussian process v (/ ). 

Modulating (14) with the Fourier modulating function <|> mfW (0 and integrating over [0,7], 
the differential equation is seen to be equivalent to the following n th order difference equation 
upon utilizing Eq. (12) of Property 2 : 

A" (im co 0 ) n Y(m)+ £ A" (im co 0 )" ~ j Y(m) = £ b } A" (im to 0 )" ~ j U(m) + e„(m). ( 1 6) 

7-1 7-1 

In this process, the equation error e(t) is transformed into the discrete frequency sequence: 

71 71+771 

e„ (m ) = A" A (im to 0 )V (m)= £ (ik (a 0 ) n ~ J c k _ m V (k ) 

0 k-m 

where V(k) is the k lh harmonic Fourier series coefficient of the white noise process v(r). 
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Interchanging the orders of summation, the preceding equation can be written as 

n+m 

e n (m)= ^a(k,m,Q a )V(k) (17) 

k^m 

where 6 a =(a ], • • a n ) and a(k ,m ,Q a ) is a frequency-dependent function of the parameters 
defined by 

n 

a(k ,m ,0 a ) = c k _ m £a, (ik co 0 )""- / , a 0 = 1 . (18) 

jM) 

As can be inferred from Property 3 , the stochastic sequence e„ (m ) is Gaussian with a corre- 
lation function that possesses finite support. The implication of this is that the covariance 
matrix for the residuals is a banded structure, specifically it is banded by the order n of the 
differential operator model, which significantly simplifies the search for the appropriate 
weighting matrix in a weighted least squares estimation. 

Defining a 2n -column parameter vector 0 as 

Q = (-a h ---a n< b h - b n y = (re a ', B b ' )' 

where prime denotes transpose, (16) can be rearranged into the following linear regression: 

Yf(»») = y(m)Q + e„(/n), m = 0, 1, ■ • M. (19) 

The 2 n -row vector of regressors y(m ) is defined by: 

y(m) = ( Y f(m), • • y %(m), y “(m), ■ ■ y £(m) ) (20) 

and the pairs ( y f(m), y J(m) ) comprising y (m), j =0, 1 n , are defined by: 

[y f(rn ), y/(rn ) j = A" (im co 0 ) n ~ J [u(m),Y(m)]. (21) 

Taking into account the fact that (19) is complex-valued, a cost function for a weighted least 
squares minimization is defined by: 

7(0) = (Y c - T C Q)'W-\Y C - T c 0) (22) 

where the following notation applies: 

Re Yo v (0) 

Re y m) 

Im y^O) 

Imy^M) 

with “Re” and “Im” meaning real and imaginary part, respectively. Assuming linearly 
independent regressors, minimizing (22) leads to the least squares estimate: 

0 = (r c 'w~ l r c r l r c 'w~ 1 y c . ( 24 ) 

This estimate cannot be biased by unknown initial or boundary conditions in the data owing 
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to the properties inherent in the modulating functions of Shinbrot type. It remains to specify 
the weighting matrix W~ l which will be done iteratively in a ‘relaxation’ algorithm that takes 
into account the parameter dependent form evident in (17) along with the white Gaussian 
assumption on the measurement noise v(/). However, the algorithm in its simplest (deter- 
ministic) least squares format: 

=(r c T c r 1 ryy c 

i.e., the estimate (24) with W=J , will first be compared with a well-known prediction error 
algorithm in order to illustrate its performance in the most elemental terms. 

3.2. A Comparison Based on Simulation 

Consider the second order system over a fixed 105 time interval: 
y(0 + 3y(r) + 8y(r) = 5 u(t), OZtZT = 105 

which possesses the transfer function: H(s) = 5/(s^+3s+8). Shown in Fig. 3 for point of 
reference is: (a) the time domain step response of the system, (b) the frequency domain mag- 
nitude plot, and (c) the effect of sampling rate on the parameter estimation errors under ideal 
zero-noise conditions. The latter is discussed below. Since there are only 3 unknown param- 
eters and n= 2, M was chosen as M=6. This means that (M +n )Fq=0.8 Hz is the highest fre- 
quency that will be extracted from the data, which is roughly the bandwidth of the system as 
seen in Fig. 3(b). This also assures adherence to the guidelines (9a) and (9b). 



TIME (s) FREQUENCY (Hz) SAMPLING FREQ (Hz) 


(a)STEP RESPONSE (b)MAGNJTUDE PLOT (c)NORMALIZED ERR0R(%) 


Fig. 3. Aspects of the simulated system 
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The output y (r) was corrupted by additive white Gaussian measurement noise. Two 
hundred Monte Carlo runs were made for each of several noise-to-signal ratios (NSR) under 
two separate conditions: (i) the initial conditions fixed at (0,0), and (ii) the initial conditions 
randomized for each run, i.e., a total of 400 Monte Carlo runs at each NSR. The input was 
u (t )=sin/ 2 /5 over the T =10s time interval for each run. All calculations were carried out in 
MATLAB. The sampling rate was fixed at 25.6 Hz, thus facilitating a standard 256 point 
DFT for the 10s time interval. 

The results are shown in Fig. 4 where “LS/MFT” means the unweighted least squares 
estimate based on the modulating function technique with the weighting matrix W chosen 
simply as the identity matrix. Also shown is the estimate based on the prediction error 
method (PEM) from the Identification Toolbox in MATLAB (Ljung, 1991). 



O 10 20 30 40 50 60 0 10 20 30 40 50 60 

(a) NORMALIZED BIAS (%) (b) NORMALIZED STD (%) 

VERSES NSR (%) VERSES NSR (%) 


Fig. 4. Output measurement noise effects for the LS/MFT and PEM. 

In the case of nonrandom (fixed) initial conditions shown in the top halves of Fig. 4, the PEM 
shows a fairly constant 6«7% normalized bias error 4 as the NSR increases from 0% to 60%, 

7 V/ 

I n e Q — 0 . * 

4 Defined as — Y\ t — • 100% where 0, is the ensemble average of the MFT esti- 

"ejill 6;* J j 

mates for the true 0y * . 
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whereas the bias for the MFT gradually increases from 0% to ~ 12% over the same range of 
NSR values. Meanwhile, the normalized standard deviations are nearly the same for this case. 
However, in the case of randomized initial conditions, see the lower halves of Fig. 4, both the 
bias and the standard deviation for the PEM show a several-fold increase and, moreover, each 
exhibits rather unpredictable behavior that is difficult to reproduce even with additional Monte 
Carlo runs. By contrast, the bias and standard deviation are very repeatable for the MFT and 
both have decreased at each NSR indicating that nonzero initial conditions have no deleterious 
effect on the MFT. 

The above comparative results have been verified on other examples as well (Fullerton, 
1991). One reason for the superior performance of the MFT is that the MFT does not have to 
estimate unknown initial conditions; another is related to the fact that the MFT is a direct 
identification technique for continuous-time models, while the PEM first estimates the parame- 
ters of a discrete-time model then converts this to a continuous-time model. Also, the low 
frequency Fourier series coefficients needed by the MFT can be computed with sufficient 
accuracy using modest sampling rates. Additional insight in this regard is gained by referring 
to Fig. 3(c) which shows the influence of sampling rate on the percent normalized enror under 
zero measurement noise conditions. Apparently the PEM needs a much higher sampling rate 
for the given 10s of data, or a much longer T time interval for the given 25.6 Hz sampling 
rate, in order to achieve a small estimation error. The influence of sampling rate on the esti- 
mation error for the MFT is much more consistent with the Nyquist sampling theorem in view 
of the frequency magnitude plot in Fig. 3(b). Finally, it is interesting to note that the MFT 
requires substantially less computer time for each set of Monte Carlo runs (about a ten to five 
fold decrease depending on whether or not the PEM has to estimate the initial conditions). 

4. Cost Functions for Aircraft Models 


4.1. Longitudinal Dynamics 

A joint cost function that reflects the equation errors for each model in (l)-(2) is 
specified as follows: 

= (t i - r^yw-HYi - iw + v(r 2 - r 2 Q 2 yw~ l (y 2 - r 2 e 2 ) (25) 


where v is a s calin g parameter introduced to accommodate the unknown variances in the two 
output measurement noises. The subscripts 1 and 2 are used to distinguish the parameters and 
data associated with the two equations (1) and (2) for this model. Specifically, 


01 = 




(26) 


where (0 a ,Q b .d b ) are the parameter vectors comprising the coefficients of 
(Aj(s), B { (s), B 2 (s )) respectively. Likewise, the regressand vectors (Y 1( Y 2 ) and regression 
matrices (r h T 2 ) are defined akin to (23) in a way that makes the vector-matrix multiplica- 
tions conformable with the defined parameter vectors. Since the parameter vector Q a associ- 
ated with A 1 {p) is coupled into each equation error, the regression matrices are partitioned 
according to 
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( 27 ) 


r i - fX fli , r^], r 2 = [r fl2 , r fc j 

such that (25) can be rewritten as 

^i(0 a .06 1 .o*p = 0^1 + r a e fl - r^e^yir 1 ^! + r a e a - r b fi b ) 

+ v(r 2 + r fl2 e a - r b e b Yw-\Y 2 + r a e fl - i\ 0*j. 


(28) 


4.2. Lateral Dynamics 

A joint cost function that reflects a weighted composite of the equation errors for (3)-(5) 
is 


y 2 (0 1 ,e 2 ,0 3 ) =£v 1 _ 1 (y, - r^jw-^Yi - r.-o,.), v 0 = i 


i-l 


(29) 


where the (v 1( v 2 ) are scaling parameters to accommodate the unknown output measurement 
noise variances and permit adjustments to the minimal values of the individual equation error 
residuals. The parameter vectors are defined analogous to (26): 



% 


i = 1,2.3 


(30) 


where (0 a ,0^ ), i =1,2,3, are the parameters comprising the coefficients in the polynomials 
A 2 (p) and (B i j(p ), B^ip )), / =1,2,3, respectively. With a partitioning of the regressor 
matrices T, akin to (27), the joint cost function (29) can be rewritten as 


'20..VW = E v i-i( y i + r « 9 « - r&yw-'or, + r^e. - r„v 

i-l 


(31) 


4.3. The AWLS/MFT Algorithm 

Equating to zero the partial derivatives of the cost functions (28) and (31) with respect to 
their arguments is a necessary condition for a weighted least squares estimate of the parame- 
ters for each model set. In the case of (28), this procedure leads to the following coupled set 
of equations for (0 a ,0^,0^): 

(r./w-ir,, + vivw'-'r^e* = - r a 'w-'Y, - vr^iy- 1 ^ 

+ r a y-'r b f> bt + viyw'-’r^ ( 32 ) 


= r b ;w-'r , + r bl 'w-'r a p a (3 3 ) 

= r b ;w~'y 1 + r„ ! 'H'- 1 r 0! e o . (34) 

Here it is important to note that each equation has been arranged such that the left hand side 
contains the symmetric normal-type coefficient matrix, premultiplying the subparameter 0 a , 
Q bi or 0 fcj , which would be used to obtain the least squares estimate of that particular 
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subparameter given the other subparameter values on the right hand side and given the 
weighting matrix W~ l . 

A basic question considered in Pearson and Shen (1993c) for n th order SISO systems is: 
How to choose the weighting matrix W1 Focusing on the partitions: 

WjR Wj 



where {1^?}^,^+.; = £[£/? +/)], etc., and the subscripts R and I denote real and 

imaginary parts, it is shown therein that W /, w =0, 1^=0, and the submatrices (W R ,W r ) are 
nearly identical, differing only in the entries corresponding to the zero (DC) frequency, cf. 
Property 3. Moreover, (W R ,Wj) are banded Toeplitz structures (the width of the band being 
n ) that depend in a known parametric way on the 0 a parameters via the frequency dependent 
function a(k ,m ,Q a ) defined in (18). Denoting this functional dependence by W=g(0), the 
exact expression of which is given in Pearson and Shen (1993c) but can be gleaned from Pro- 
perty 3, the SISO “adaptive weighted least squares” algorithm iteratively seeks a solution for 
the pair (Q^wzj*^) from the equations, cf. (24): 

0AWLS = (ryw'- 1 r c r 1 ryw'- 1 y c (35a) 


W=S(W>)- (35b) 

The initial guess for starting the iterations is either W=I , or the weighting matrix that arises 
in the ideal situation in which the equation error e (t ) is just a Gaussian white noise process, 
i.e., not dependent on the 0 parameter. The resulting solution to (35a) is inserted into (35b) 
thereby obtaining a new estimate of W which is reinserted into (35a), etc. This “relaxation” 
type algorithm has been used by many researchers as a tool to achieving a weighted least 
squares estimate that is close to a maximum likelihood estimate under certain conditions. It 
almost always converges though there is no known proof of this factorage. 

The following point should be emphasized in comparing the various estimates available 
under the MFT framework: experience shows that the AWLS estimate is far less sensitive to 
the chosen bandwidth F B than the ordinary (unweighted) least squares estimate. As evidence 
of this statement, a fourth order SISO model was used to relate the longitudinal pilot stick 
command to the body axis pitch rate q(t), the transfer function of which had a bandwidth 
around 0.5Hz. During the modeling process, a set of modulating function bandwidths ranging 
from 0.3Hz ~ 1.0Hz (A/ ranging from 12 to 40) was chosen for the LS/MFT and AWLS/MFT 
algorithms. An in-between type algorithm WLS/MFT was also used for comparison; this is 
the algorithm obtained by fixing the weighting matrix W to that which corresponds to the 
ideal white-noise residuals case mentioned above. In this case, the weighting can be obtained 
in closed-form fashion from Property 3 by substituting P (s )=1 and utilizing the identity: 

_ (-l)'(2n)! 

t ir*/ ' («-/)!(»+/)! • 

The results are summarized in Table 1 where the goodness-of-fit ratio, s/e, is defined in (36) 
below. The ‘constrained’ and ‘unconstrained’ columns under AWLS/MFT refer to use of this 
algorithm in a mode whereby the estimated model is forced to have stable poles (constrained), 
or not (unconstrained). This constrained option for AWLS is invoked by checking the poles 
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of the estimated model during the iterative process and shifting any right half plane poles to 
the left half plane, via reflection about the imaginary axis, before continuing the iterations in 
case that event occurs. This can occur for aircraft modeling since the systems are frequently 
lightly damped. Concerning the results in Table 1, it is seen that the models for the LS/MFT 
and WLS/MFT are either stable with poor s/e ratios, or unstable with unacceptable s/e’s. By 
contrast, the estimated systems for both constrained and unconstrained AWLS/MFT 
algorithms yield much more consistent time domain performance with no drastic variations 
like the LS/MFT and WLS/MFT algorithms, regardless whether the estimated models are 
stable or not. Hence, these results help establish the claim that the AWLS/MFT algorithm 
not only performs better, but is less sensitive to the pre-chosen modulating bandwidth F B . 


algorithm 

F b (Hz)/M 

AWLS 

constrained 

AWLS 

unconstrained 

WLS 

LS 

0.3/12 

stable 

s/e = 5.84dB 

stable 

s/e = 5.84dB 

unstable 
s/e = - 4.78dB 

stable 

s/e = - 6.79dB 

0.4/16 

stable 

s/e = 6.83dB 

stable 

s/e = 6.83dB 

unstable 
s/e = 2 09dB 

unstable 
s/e = - 3.39dB 

0.5/20 

stable 

s/e = 7.13dB 

stable 

s/e=7.13dB 

stable 

s/e= 1.86dB 

unstable 
s/e = 1.53dB 

0.6/24 

stable 

s/e = 6.91dB 

unstable 
s/e « 7.56dB 

stable 

s/e = 3.08dB 

unstable 
s/e = - 7.41dB 

0.7/28 

stable 

s/e = 6.51dB 

unstable 
s/e = 8.07dB 

stable 

s/e = 3.32dB 

unstable 
s/e = - 53.22dB 

0.8/32 

stable 

s/e = 6.41dB 

unstable 
s/e = 8.10dB 

unstable 
s/e = - 27.63dB 

unstable 
s/e = - 41.04dB 

0.9/36 

stable 

s/e = 6.56dB 

unstable 
s/e = 7.94dB 

unstable 
s/e = - 33.98dB 

unstable 
s/e = - 27.60dB 

1.0/40 

stable 

s/e = 6.47dB 

unstable 
s/e = 8.00dB 

unstable 
s/e = - 61.97dB 

unstable 
s/e = - 121. ldB 


Using physical flight data to build a fourth order model linking the longitudinal pilot 
stick movement (input) and the body pitch rate (output). 


Table 1 . Sensitivity of MFT algorithms to the chosen modulating bandwith 



























































The following relaxation type algorithm is an extension of the algorithm just described 
for solving (35) but aimed at solving the set of equations (32) - (34), together with the equa- 
tion for W that corresponds to (35b) in this case, over the parameter set (O^O^.O^.W): 

The AWLS/MFT Algorithm for the Longitudinal Dynamics Model (l)-(2) 

1. Pick a scaling parameter value for v, v^O, and estimate an initial value for the pair 
(0 a » W ) through the model: 

A&XaQ) + <?(/)] = [B 1 (p)+B 2 (p)]u(t) 

using the SISO AWLS/MFT algorithm relative to (35). 

2. Substitute the values for (0 a , W) from step 1 into (33) and (34), then solve for the pair 

(0*.. 

3. Estimate a new 0 a from (32) using the values for (0 fc| , 0^, W) from the previous step. 

4. Check if the parameter value for the new 0 a has changed or not, based on a user-chosen 
percent change in norm. If yes, compute a new value for W from (35b) and go back to 
step 2, otherwise stop. 

5. Check the system output-signal-to-output error ratios S/E, defined below in (36), for each 
of the two models (1) and (2) to see if they are in rough agreement with each other. If 
not, try a new value for v and repeat steps 1 ~ 5. 

The algorithm for the lateral dynamics is similar to the above, i.e., the partial derivatives 
of (31) with respect to its arguments are equated to zero and arranged in a manner analogous 
to (32) - (34). Further details are available in Shen (1993), including a computationally 
efficient algorithm for inverting W (total flops of order M 3 ) that exploits its banded structure 
at each step and is very robust for large ( M ,n ). 

4.4. Modeling Results 

The results of applying the algorithm to flight data for a high performance aircraft will 
be discussed in this section and, whenever possible, comparison will be made with the model- 
ing results from an established maximum likelihood/output error algorithm which is based in 
the time domain. (See description of the latter algorithm in Klein (1993).) The figure of merit 
used to assess the quality of the resulting model in relation to the given data is the output- 
signal-to-output-error ratio, S/E, (in decibels) as defined by 

S/E = 20 iogio ^f^j , e{t) = y(t) - y(t), 0<J<T (36) 

where RMS(y ) means the root-mean-square value of the output signal y{t) over the [0,T] 
time interval, and y (t ) is the estimated output using the model. 5 Subscripts on S/E are used to 
distinguish a particular output among the multivariable channels. 


5 When required, an observable state equation model is employed and the unknown initial condi- 
tion for generating y(t) is estimated from a Luenberger observer r unnin g backwards in time driven by 
the given input/output data over the [0,T] time interval. 


- 16 - 



The data for the longitudinal case and the corresponding modeling results are given in 
Fig. 5 and Table 2, respectively. 






Fig. 5. Input/output data for the longitudinal dynamics, and 
superimposed model output responses. 

The parameter symbols in Table 2 are defined by the transfer function relations: 
a(s) _ ^ a 5 + £ a q(s) _ A? 5 + Eg 

“CO s 2 + a:^ + a: 0 ’ “CO ~ s 2 + a:^ + k 0 

where u can be either 6/, or depending on whether it is the transfer function obtained by 
using as input the signal “inside” the loop, referred to as ‘Open Loop’ in Table 2, or “out- 
side” the loop, referred to as ‘Closed Loop’ in Table 2. (See Fig. 1.) 


Parameter 

Open Loop 

Closed Loop 


ML 

AWLS/MFT 

AWLS/MFT 

KO 

0.695 

0.698 

0.494 

Kl 

0.360 

0.324 

0.995 

Aa 

0.002 

-0.065 

-0.013 

Ba 

-1.558 

-1.523 

0.162 

Aq 

-1.479 

-1.488 

0.143 

Bq 

-0.328 

-0.273 

0.043 

(S/E) a (dB) 




(S/E) q (dB) 

| 

1 

| 


Table 2. Modeling results for the longitudinal dynamics. 
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The symbol ML in Table 2 refers to the maximum likelihood/output error identifier described 
in Klein (1993). Clearly, the S/E ratios for the Open Loop case are nearly the same for each 
method in each of the respective output channels. The same is true for the Closed Loop case, 
although only the AWLS/MFT values are listed for the sake of brevity. Moreover, the meas- 
ured and predicted model responses are quite close to one another, as evidenced by the right 
hand sides of Fig. 5. 

The scaling parameter value used for the AWLS/MFT results in Table 2 was v=0. The 
results of using other v values for the Closed Loop case are depicted in Fig. 6 showing the 
tradeoffs that can be obtained by weighting the two output channels differently. The ‘itera- 
tions’ referred to in Fig. 6 count the number of loops passing through the AWLS/MFT algo- 
rithm listed in the previous section. Ordinarily one should expect the curves of the type 
shown in Fig. 6 to be monotonic, leveling off to some asymptotic values for each channel. 
This has been the case for other studies (Shen, 1993), but here it possibly indicates a straining 
of the modeling assumptions, e.g., linearity. 



Fig. 6. Effect of different v values on S/E ratios for the longitudinal dynamics. 

The data and resulting transfer functions for the lateral dynamics study are shown in Fig. 
7 and Table 3, respectively. Only the AWLS/MFT technique was successful in yielding a 
useful result in this case, i.e., the ML technique failed. The scaling parameters were both 
chosen as unity for the AWLS/MFT, i.e., vpv^l. As seen at the bottom of Table 3, reason- 
able parity was achieved in the S/E ratios for each of the three output channels; hence, there 
is no compelling need to change the v values from unity in this case. Again, the measured 
and predicted model transient responses are in good agreement with each other as shown by 
the three plots on the right hand side of Fig. 7. 
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Fig. 7. Input/output data for the lateral dynamics, and superimposed model output responses. 

Table 3. Modeling results for the lateral dynamics. 

P _ 1.8131s 2 +.3099s+.2001 p .0082s 3 +.3701s 2 +.0819s+.0543 

S a s 4 +. 31 08s 3 + 3. 2640s 2 +. 2059s+. 2536 S t A 

p _ 3.0633s 3 -2.2515s 2 -1.2644s-.0049 p _ .1909s 3 -.S163s 2 -1.2768s+.1098 
$a A 

r -.1627s 3 -.1845s 2 -1.3654s-.0148 

s a A S r A 

(S/E)p = 9.12 (dB) (S / E) p = 12.31 (dB) (S/E) r = 9.96 (dB) 


A 

-.2031s 3 -. 1630s 2 -.791 2s-. 0962 
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5. MFT Extensions 


5.1. Frequency Transfer Function Analysis 

Methods for determining the transfer function in the frequency domain for a stable linear 
system from input/output data include classical correlation and spectral analyses, as well as 
the direct ratio of Fourier transforms and steady state sinusoidal measurements. Each of these 
“nonparametric” identification techniques require either a statistical stationarity assumption 
on the data, or a periodic steady state condition to be established, before initiating calculations 
of the transfer function at pertinent frequencies (Ljung, 1987). Notwithstanding noise con- 
siderations, long data lengths may be required in order to achieve good accuracy due to the 
stationarity or steady state assumption, or to eliminate end point effects in case a direct ratio 
of Fourier transforms is used on time-limited data. By way of contrast, the modulating func- 
tions (10) will be used in this Section to extract the frequency content in short data lengths in 
order to set up a least squares estimation of the transfer function at selected frequencies. 
Since short data lengths are used there is no assumption of steady state operation or stationar- 
ity of the data, though there must be present sufficient energy content in the data at the 
specified frequencies in order to avoid degeneracy in the least squares estimate. 

Consider a stable linear system with two inputs (u(t),v(t)) and a single output y(t) 
modeled by the equation: 

A(p)y(t) = B(p)u(t) + C(p)v(t) + e(t) (37) 


where (A(p)Ji(p),C(p)) are polynomials in the differential operator p=d/dt of degree less 
than or equal to an a priori integer n, and e(t) again represents the effect of modeling errors. 
Given the input-output data [n(r),v(r),y(f)] over a finite set of time intervals 
{[t r ,t r +T], r= 1, • • N}, the problem considered here is to estimate the transfer functions 
G (/ (a)=B (z co )jA O' co) and H 0 co)=C (/ c 6)1 A (/' co) at a finite set of frequencies 
{ 7 co 0 , 7=0,1 ■ • M+n }, 6 where coq is the user selected ‘resolving’ frequency and M a chosen 
nonnegative integer. The time intervals are each chosen of length T =2n/co 0 and need not 
necessarily be disjoint. 

As in the previous formulations, the model (37) is rearranged into the equation error for- 
mat and modulated with ^ mjn (/). Thus, with respect to the time interval [t r ,t r +T]: 

T 

1(0 

Invoking Property 2 , Eq. (12a): 

m+n r 1 m+n 

£ c k-m \ A ( ik <o 0 )y* [r ] - B O'* <o 0 )U k [r] — C (ik co^V* [r ] = £ E [r ] (3g) 

k~m ^ k-m 

where the following notation applies to the Fourier series coefficients in the above (Z can be 
U ,V ,Y or £): 


A(p)y(t+t r ) -B(p)u(t+t r ) - C(p)v(t+t r ) dt = ^ mjl {t)e{t+t r )dt. 


6 Thus, assuming the system bandwidth (£>g is known, the choice in (M,co 0 ) such that 
(M+n )co0~co 5 covers the bandwidth at the knots k cdq, 7=0,1 • • M+n. 
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1 T 

Z*[r] = —jz(t+t r )coska)Qtdt - i— |z (/ +t r )sink (atfdt 


Z£[r] 


- i Zftr]. 


( 39 ) 


Define parameters (a k R ,a. k f ,fi k R fid ,y k R >Yk) f° r the real and imaginary parts of the polyno- 
mials {A {ik co 0 ),B (ik a> 0 ),C (ik co 0 )), each multiplied by c k _ m , e.g., 

c k-rn A = o- k R + ia k l . 

Values of the transfer functions G (i to) and H (i to) are related to these parameters: 


G (ik (o Q ) = 


+ i$k 
a k R + i a k 


H (ik coq) - 


y** + nj_ 

a k R + i a k 


The DC value can be handled by normalizing A (0)=1 so G (Q)=B (0) and H(0)=C(0). With 
this normalization and the notation of (39), the real and imaginary parts of (38) can be rear- 
ranged into the following standard linear least squares equation format represented by stages 
according to the values assigned to m . 

Initial Stage (m=0): 


Y 0 [r] = DVttrlQ* + £<)[r], r- 1,2 ■ • N 

km 0 


(40) 


where the data-related quantities are defined by the vector-matrix equations: 

[t/&[r] V c Q [r] 

I nr-IVl — 

0 


Y 0 [r] = 


Vo[ r ] = 


0 


0 


= 


Y k [r ] Y{[r] U c k [r ] £7/[r] V c k [r] Vftr] 

-Y s k [r] Y k [r] -Uftr] U c k [r] -Vftr] Vftr] 


(k*l) 


and the Q k parameters are defined by 

\g( 0) 


0n = 


H( 0 ) 


, 0* = col (~a k , -a/, pf, p/, y k , y/ ) (££1). 


Stage m (m=l,2 • M): 

Again, a 2 dimensional vector equation is derived from the real and imaginary parts of 
(38) for the case m^l with the same notation as above: 

Y m [r] = v n+ m ['•]9n+m + ^ I'l r= 1,2 • • N (41) 

where Y m [r ] is defined in terms of the parameters and the data of the preceding stages by 

n+m-\ 

Y m [r] = - £ vd r Wk- 

kmm 
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The residuals for all stages are given by 


n+m 

1 — Sr c k-m 

k*m 


mr] 
-mr] ' 


Remarks, (i) Comparing (40) and (41), it is seen that the m=0 (initial) stage carries the 
major computational burden since there are many more parameters to be determined at this 
stage, i.e., (2+6/i) for (40) verses just 6 for each succeeding stage in (41). Notice that the 
same data can be used for each stage and that N must satisfy: 2 N >(2+6 n ) in order to have 
more equations than unknowns, (ii) In a general multivariable situation, the transfer function 
representation for each output yj, 7 =1,2 • • m y : 

™ u m u B it (s ) 

yj(s) = t,Hjk(s)u k (s) = Z-r—7-ru k (s) 

k- 1 i»l A jkV) 

facilitates the following model for each j : 

77l„ 

Aj(p)yj(t)= i,B jk (p)u k (t) + e(t) (42) 

k-1 


m u 

where Aj (s jk ( s ) and B jk (s) is defined as the polynomial: 
k = 1 

BAs) 

Bjk(s) = T^Ts ) Aj(s) - 

Therefore, the multi-input single output formulation developed above can be applied to each 
output equation (42). The fact that the pairs (Aj(s), Bj k (s)) t k=l, • • m u , are not generally 
coprime seems not to be a problem under low measurement noise conditions since it is the 
ratios B jk (i aS)!Aj (1 co) = B jk (i w)IAj k (i 00) that are sought at the specified frequencies. This 
avoids a difficult issue of sufficient parametrization and minimality for a state space model 
which is a separate issue. 


5.1.1. A MIMO Example 

As an example, consider the system with the transfer function relations (notice that 
H 2 (s) is high-pass): 

y^s) = G^sMs) + G 2 Cs)v(s) 

_ 12* 2 +4875+582 + 0.7j 2 +157.5s+504 

~ s 3 +65s 2 +456s+ 1978 1 5 3 +9s 2 +455s+881 ^ 


y 2 (s ) = H i(s ) u ( s ) + H 2 (S )v (s ) 
2s+160 


s 2 +205+160 


. . 5 2 +50s+54 , . 

u(s) + — v (5 ). 

s 2 +40s+500 


The system was simulated over a 305 time interval using as inputs: 

u (/) = a random binary signal 
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V (r ) = sin(25/ +0.9434) + sin(6f) + sin(20r 2 ) + sin(5r 3 ). 


Order n=6 modulating functions were utilized to accommodate the degrees of the poly- 
nomials {A )) i* 1 (42) for the above transfer functions. Choosing the desired 

bandwidth to be covered and at what resolving frequency led to the goal of estimating the 
four frequency functions at the knots { k co 0 , 7=0,1 -15} where coo=27t/4=1.57 rad/sec, i.e., 
(o B =23.55 rad/sec and T time intervals of 45 duration. Taking into account the number of 
unknowns for the initial stage, i.e., (2+6n)=38, twenty-seven [t r ,t r + 4] time intervals were 
chosen with 1 s overlaps, i.e., [t r ,t r +T] = [r,r+4], r=0,l • • 26, which facilitates about 50% 
more algebraic equations as unknowns upon which to base the least squares estimate for the 
initial m=Q stage. The sampling frequency was 200 Hz in utilizing the DFT to calculate the 
Fourier series coefficients of the data on each 45 subinterval. In the absence of noise, the 
estimates of the frequency functions Gj(i co) and Hj{i co), 7=1,2, are essentially perfect at the 
knots {7 co 0 , 7=0,1 • • 15}, co 0 =2n/4=1.57rad/sec, as expected. In order to test the sensitivity 
to noise, 50 Monte Carlo runs were made at each of several noise levels with white Gaussian 
additive noise at the two outputs. The ensemble mean magnitude and phase plots are shown 
in Fig. 8 for G^/to) and G 2 0’co), and in Fig. 9 for //jO'ca) and H 2 (i to), with magnitude on 
the top half and phase on the bottom half for each transfer function. Also shown for each 
transfer function is: (i) a low noise condition (SNR = 35 db) shown on the left hand side of 
each figure, and (ii) a higher noise condition (SNR = 5 db) shown on the right hand side of 
each figure. Generally, the estimates are better at the low frequency end of the spectrum as 
might be expected for white measurement noises. 


5.2. Linear Time Varying Differential System Models 

Consider the linear time-varying differential operator model 

P n y{t) + 'Za j (t)p n ~jy(t ) = Y,bj{t)p nu ~ j u(t) + e(t) (43) 

i m i 7-1 

where the variable coefficients are parametrized by linear combinations of smooth functions. 
To facilitate generality, it is convenient to adopt the row and column vector notations, “<” 
and “>”, such that each coefficient is represented by an inner product: 

o.j>, b n _j(t) = < gj (t) $j> ( 44 ) 

where <fj (t ) and <gj(t ) are given row vector-valued smooth functions, and (a ; >, (3 j>) are 
column vector-valued parameters. Stacking the latter columns into a pair of column vector- 
valued parameters (co, |3>), viz. 

oo = col (ao>, ai> ■ ■ a n _i>), p> = col (Pq>, • ■ P„ u _i>) 


and utilizing the Leibniz formula in order to achieve a differential operator model amenable to 
Property 2, it is found that (43) together with (44) is equivalent to the model: 



n - 1 . 


n u - 1 

p n y(t) + 

Zp < WO 

OO = 

£/>*<«*( r) 




*«o 


P> + e(t) 


(45) 


where and <y k ( t) are row vector functions related to time-modulated versions of the 

input/output data, modulated in the time domain by derivatives of the given functions <f j(t) 
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Ptiaae Responae(degree) 


(a) G i(i to) 


Mem Values for SNR=35(dB) Case 


Mew Values for SNR=5(dB) Case 




Mew Values for SNR*35(dB) Case Mew Values for SNR*5(dB) Case 




(b) G 2 (/co) 


Mean Values fa SNR=35(dB) Case Mean Values for SNR-5(dB> Case 




Mean Values for SNR*35(dB) Case 


Mean Values far SNR*05(dB) Case 



Fig. 8. Measurement noise effects for Gj and G 2 
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(a) H i(i co) 


Mean Values for SNR*35(dB) Case 


Mean Values for SKR=5(dB) Case 


(b) H 2 {i CD) 




Mean Values for SNR=35(dB) case 



Meao Values for SKR«5(dB) case 



Mean Values for SNR=35(dB) Case 



Mean Values for SNR=35(dB) Case 



Mean Values for SNR=5(dB) Case 



Mean Values for $NR*5(dB) Case 



Fig. 9. Measurement noise effects for H y and H 2 
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and in accordance with the following definitions: 

<u k (t) = row(<v 0k (t ),<v lk (t ) • • <v nu _ l k (/)) 

<y k (t) = row(<z 0i .(r),<z u (r) • • <z„_ u (0) 
and 


<v jk ( t) = u(t )P jk {p )<gj (r ), <z jk (r) =y(t)P jk (p )<fj (t ) 
with Pj k {p ) defined as the differential operator: 


Pjk(P) = \ 


0 for j< k 
(“IV* \{\p j ~ k for £ 


Applying Property 2 , Eq. (12), to (45) leads directly to the regression equation (19) where 
y(m )0 = A" 


L (im co 0 )* <Y k (m ), £ (im (£> 0 ) k <U k (m ) 


k ~ o 


*«0 



-a> 


.p> . 


5.2.1. A Variable Damping Example 

The following second order system with a parabolically varying damping term and 5 
unknown parameters was simulated under a variety of conditions: 

y(t) + (£ +c 1 /+c 2 t 2 )y(0 + ay(t) = bu(t), 0<,t<T=l0s. 

Via the Leibniz formula, this model is equivalent to 

P 2 y ( 0 + C py(0 + Cl [p(D’CO) -j(o] + c 2 \p(t 2 y(t)) -2/y(r)j + ay(t ) = bu(t ) 

which is an illustration of (45). The input was u (r)=sinr 2 /3 for each run, and the parameters 
(a ,b ,c \,c 2) were fixed at the values: (8, 5, -1, 0.1), respectively. Two cases were studied for 
the £ parameter: £=3 in which the “frozen poles” of the system are stable for each re[0,10s], 
and £=2 in which the “frozen poles” are unstable for approximately 50% of the 10s time 
interval. (The parabolas for these two cases are shown in the top halves of Figs. 11(a) and 
(b) respectively.) Each £ case included: (i) 300 Monte Carlo runs at each of several noise- to- 
signal ratios with white additive measurement noise on the output signal y (r), and (ii) struc- 
tural errors in the modeling of the damping term by estimating first a linearly time-varying 
damping term, and then a constant damping term. The data was noise-free for this aspect of 
the study. 

The noise effect results are summarized in Fig. 10 where the left side shows the normal- 
ized bias errors and the right side the standard deviations. The results are more or less the 
same for the two £ cases and demonstrate that the estimation technique has similarly good 
noise-immunity properties as exists for time-invariant systems. 

The effect of structural errors is shown in Fig. 11 where the left and right hand sides 
pertain to the t=3> and £=2 cases, respectively. In the absence of noise, the parameter esti- 
mates are essentially perfect, as expected, and yield the parabolas for the two cases as shown 
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in the top halves of Fig. 11. With structural modeling errors, the estimates of the damping 
terms are as indicated in the top halves where it is noted that modeling the damping as a con- 
stant resulted in an unstable system for the £=2 case. Using a different input signal : 
u(0=sinf 4 /4, i e., different from that which was used to estimate the parameters, the resulting 
outputs of the estimated models are compared with one another in the bottom halves of Fig. 
11. As expected, the more severe the error in structural modeling, the greater the error in the 
ability of the model to predict the response, though the differences are perhaps not as large as 
might be anticipated. 



(a) NORMALIZED BIAS (S) (b) NORMALIZED STD (2) 

VERSES NSR (%) VERSES NSR (7.) 


DAWPING VS TIME (») 




ESTIMATED MODEL RESPONSE VS TIME (s) 



Fig. 10. Measurement noise effects. 


Fig. 11. Structural error effects. 


5.3. A Class of Nonlinear Input/Output Models 

Parametrized state vector equations of the normal form: 

x =f(x,u,d) 

y =h(x,u ,6) (46) 

constitute the most common starting point for methods aimed towards the parameter 
identification of deterministic continuous-time nonlinear differential systems. Given the 
input/output data [n (r ),y (t )] on some time interval, (Kt <J , together with the parametrized 
functions f(x,u,Q ) and h(x,u,Q), the known methods include: quasilinearization, invariant 
embedding, state variable filters, and a variety of nonlinear filtering methods such as stochas- 
tic approximation and extended Kalman filtering. These methods are invariably implicit in 
terms of the cost function dependence on the parameter vector 0. By contrast, if an 
equivalent parametrized input/output differential operator equation exists in the following 
equation-error format (p=d/dt): 
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( 47 ) 


£ L gy (&)Fjk (“ O’ )Pjk (p)E k (u,y) = 0, g 0 = 1 

y-o *«i 

then it will be shown below that given the input/output data over [0,7] it is possible to 
specify an explicit cost function of the form 

J 0) = £ E r jk 8j (6)S* (6) (48) 

y- o *-o 

possessing the properties: (i) J (0)^0, with 7 (0)=O for any value of 0 that satisfies the model 
(47), and (ii) the rj k can be computed via integral operations on transient data over the obser- 
vation time interval [0,7] without the need to estimate unknown initial conditions. Hence, 
minimizing J (0) in order to ameliorate modeling errors leads to a one-shot least squares esti- 
mate of the system parameters, which may or may not be unique depending on the adherence 
of appropriate nondegeneracy conditions on the data. 

The formulation leading to (48) is a continuation of earlier developments in Pearson 
(1988, 1989) relating to the least squares parameter estimation for input/output models like 
(47) based on Shinbrot’s classical method of moment functionals. The formulation in this 
Section is taken from Pearson (1992) and exploits the Property 4 delineated earlier. 


5.3.1. Modeling Considerations 

As a first step it is assumed that the system model can be arranged into the form (47) 
where the gj (0) are given functions of the parameters 0, the Pj k (p) are fixed polynomials of 
degree n in the differential operator p=d/dt, and the (E k (u ,y)J r j k (u ,y)) are specified func- 
tions of the pair (u ,y). 7 Starting from a state equation model (46), it is certainly not the case 
that all such models admit to an “external” differential representation (Nijmeijer and Van der 
Schaft, 1990) much less the special form (47). A simple example illustrating this model is the 
differential equation for a unit mass particle subject to a given force u (t ) and with drag pro- 
portional to the velocity squared: 

y +0i(y) 2 -0 2 w =0 (49a) 


where y (t ) is the displacement of the particle. Utilizing the differential identity: 
p 2 (y 2 )=2yp 2 y +2(py ) 2 , (49a) is equivalent to the following differential operator equation 
which is of the form (47): 


p 2 y +Q i ( i Ap 2 y 2 -yp 2 y) -02 w = 0- 


N _ . 

Higher order identities can be employed if a power series like £0y(yV, N>2, 

j-i 


model the dissipation term in (49a). 


(49b) 
is used to 


In some cases we can obtain a model where all the Fj k ’s in (47) are constants, i.e., the 
more restricted form: 


7 Sufficient smoothness and stability properties are tacitly assumed to assure the existence and 
uniqueness of bounded solutions y(t) satisfying (47) given any admissible 0, initial conditions and 
bounded inputs u (t) on [0,7]. 
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( 50 ) 


E Y,gj(Q)Pjk (p ) E k ( w .v ) = o, g 0 = !• 

jrn 0 *-1 

In such cases the model is said to be exact, a term which is suggested by the notion of an 
“exact” or “total” differential in the calculus because the equation error for such models can 
be integrated exactly given the data on [0,r], An example of this class is the following two 
compartment chemical kinetics model from Walter (1982): 

x 1 = -0^2 -02(1-03*2)* i +u 


*2 = 0 2 (1 -03* 2> x l -04*2 

y =*i- (51a) 

Assuming the output y (t )>0 for all t, the nonmeasured state *2 can be eliminated from the 
above model thereby arriving at the following equivalent input/output representation: 

/> 2 (lny) -p{—) +0 1 0 2 03>’ + 0 2 03 (py ~ u ) +04(p(l n 3') “ — ) +(0i +0 2 )04 = 0- (51b) 

y y 

As an important modeling consideration, which applies independently of the method used 
for identification, it was shown in Pearson (1989) that a necessary condition for a well -posed 
problem involving (47) is that the vector function (go(0)*£i(0) ' ' £ n ,(0)) be an injective, i.e., 
single-valued, function of the parameters else the parameter identification problem will be 
plagued by nonuniqueness. Presuming this condition is upheld and go is normalized to unity, 
it is, therefore, assumed that the function g (0) defined by: g (6)=(g k (Q), ■ ■ g Bl (0)) is injective. 


5.3.2. The Cost Functions for Least Squares Minimization 

Consider first the class of exact differential models (50). Modulating (50) with <J> m „(r), 
integrating over [0,r] and applying Property 2 , we obtain the modulated equation error 
function e (Q,m ) defined by 


e (Q,m ) = J^gj (0)yy (m ) (52) 

y'-o 

where the y y ( m ) play the role of regressors and are defined by 

"2 

Y ) (m)~ A" £ Pjic ( im a> 0 )V k (m ) (53) 

*-l 

and the V k (m) are defined by 

T 

V k (m) = jf^(u(r),y(r))e" imt0o, dr. (54) 


Forming the norm of e(6,m) over the frequency range [0,M], the cost function for a one-shot 
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least squares minimization is defined by (* denoting complex conjugate): 8 

n i M 

J(0) = L r jk = £Yy( m )Y**("0- (55) 

it«o /n=o 

Consider next the more general class of inexact input/output models represented by (47). 
Modulating (47) with <|) m>n (0. integrating over [O.T] and employing Property 4 , we obtain 
the following expression for the equation error relative to this class of inexact differential 
operator models: 

? (0,m ) = 'Lgj (9)7 j ( m ) (56) 

> 0 

where 

Y/( m) = £ W jk( m y®yjk( m ) (57) 

*-l 


Y;* On) - A n P jk (im a Q )V k (m ) 
and the (V k ( m ),W jk (m )) are Fourier series coefficients of data-related functions defined by 

i r 


V t (m) = jr£ t (u(l),y(0)i-"'" 


1 


W li (.m) = jtF jl: (u(t),y(t»e 


-/moot/ 




(58) 

r 

(59) 

(60) 


The least squares cost function for a one-shot estimate in this case is defined by, c.f. (55), 

•f(0) = £ £g, (.0)8 k (0)r jk , fj k = £ yj ( m )y k *(m). (61) 

jm 0 k**Q m* 0 

To illustrate the above notation for the inexact model (49b), the equation error for this 
linear regression problem is given by 

? (0,m ) = y 0 (m ) + ©jY^m) + 0 tf 2 ( m ) 
where the regressand sequence Yo (m ) is defined by 

Yo (m ) = A \im co 0 ) 2 Y (m ) (62) 

and the two regressor sequences (Yi( m )»Y2( m )) by 

Y \(m) = ViA 2 (im co Q ) 2 Y 2 (m ) - Y(m y»A 2 (im co 0 ) 2 y (m ), y 2 (m ) = -A 2 U(m) (63) 

where Y 2 0n) is defined by 


8 Noting that M in (55) actually corresponds to the frequency M C0q suggests that in choosing the 
pair (M ,co 0 ), their product should correspond roughly to the system bandwidth. 
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Y 2 (m) = 



im copf 


dt. 


5.3.3. Implementation and Examples 

The flow of calculations involved in setting up a least squares estimate for implementing 
the above equations is shown in Fig. 12. Assuming an accurate computation of the data- 
related functions E k (u(t),y(t)) and Fj k (u (t ),y(t)), the Fourier series coefficients in (59) and 
(60) can be carried out efficiently using DFT/FFT techniques as discussed in Section 2.1. 
Notice that for a finite M , only a finite number of such coefficients need be computed in the 
case of exact differential operator models, but a larger number is needed for inexact models in 
order to implement the linear convolutions in (57), the actual number of which depends on the 
degree of smoothness in the functions (E k (u (t),y(t)), Ej k (u (t),y(t ))) on [0,7* ] . 



Fig. 12. An explicit least squares estimation for nonlinear systems 
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Example 1. A simulation study was carried out for the inexact model (49) using data gen- 
erated by the input signal w(r)=sinf 2 /5 over a T =105 time interval with zero initial conditions. 
Thus, the resolving frequency was 1/7" =0.1 Hz. The sampling rate was 25.6 Hz, and 
Simpson’s rule (together with a 256 point DFT) was used to calculate the Fourier series 
coefficients of the data-related functions: (u (r),v (0.y 2 (0). O^f^lO. Zero mean white Gaus- 
sian noise of varying degrees of intensity was added to the output signal with 50 Monte Carlo 
runs made at each noise level. Since the parameters enter into the model linearly, the least 
squares problem is a linear regression once the regressand and regressor sequences are calcu- 
lated from (62)-(63). With only 2 unknown parameters and noting that these sequences yield 
2 values for each m, i.e., the real and imaginary parts, the sequences in (62)-(63) were com- 
puted for m =0,1,2 yielding 6 real equations upon which to base the least squares estimate 
corresponding to each input/output data pair. The mean and standard deviations of the param- 
eter estimates for the 50 Monte Carlo runs at each of 5 noise levels are plotted in Fig. 13. 
These results demonstrate that the estimator gave values quite close to the true parameters 
on-average, but the variance for the dissipative coefficient 0j was significantly larger than that 
of the 0 2 parameter. 

Figures 14 and 15 contain some additional information relating to the simulations of 
Example 1. Figure 14 shows output data for one run at the 12.8% noise level together with 
the simulated model output using the mean values of the ( 0 i, 02 ) parameters estimated at this 
noise level for the 50 Monte Carlo runs. This shows that on-average the bias in the output 
signal is quite small in relation to the system output signal, which is to be expected in view of 
Fig. 13. Figure 15 is meant to indicate the ability of the model to predict the system output 
under ideal no-noise conditions, using the parameters estimated from the noisy data at the 
various noise levels. The figure-of -merit used to illustrate this is the same as defined in (36). 

Example 2. A second simulation study was carried out for the exact model (51). Some con- 
ditions were the same as for Example 1, i.e., u (t )=sinr 2 /5, O^r^lOs, but the initial conditions 
were taken as y(0)=y(0)=2, which resulted in output data well above the unity value 
corresponding to which difficulties with the functions lny (t ) and u(t)/y(t ) could be expected 
in view of the model (51b). Although the parameters enter nonlinearly into this model, there 
is basically a one-to-one map between the 0 and g (0)=9 parameters so that this problem is 
also a linear regression in the 9 parameters. As in Example 1, zero mean white Gaussian 
noise of varying degrees of intensity was added to the output signal with 50 Monte Carlo runs 
made at each noise level. The mean and standard deviations of the estimated values are tabu- 
lated in Tables 4 and 5 for the 9 and 0 parameters, respectively. The nonlinear mapping from 
9 back into the 0 parameters accounts for the greater variability in estimating the model 
parameters. In spite of this greater variability, the average of the noise-free model responses 
using these parameter estimates exhibit reasonably good performance as shown in Fig. 16. 
Thus, if 0(/), /'=1,2 • • 50, represent the estimated 0 parameter vectors obtained from the 50 
Monte Carlo runs at a particular noise level, and y(t,Q) the noise-free model output 
corresponding to any particular 0, then a time response in this figure was obtained from the 

ensemble average: l/50£y (r ,0(/ )). A measure of the variance in these time responses is 

i - 1 

summarized in Fig. 17 which gives the standard deviations at each instant of time of the 
model responses y(f,0(/)) over the ensemble for i =1,2 • • 50. 
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Example 2 Simulations 




6,0, 

e. 

e,(e,+e,) 



true values 

0.246 

1.170 

2.000 

3.02 

noise level 

S/E (dB) 


0.247 

1.164 

1.974 

2.988 

■igiBiiia 

46.97 

ftd 

0.0776 

mmm 

0.2477 

0.1951 

0.16% 

9.100 

mean 

0.227 

1.115 

1.941 

2.925 

IM 

29.29 

std 

0.3734 

0.5309 

1.1869 

0.9314 

0.79% 

11.076 

mean 

0.193 

0.918 

1.509 

2.335 

-rgimaoi 

19.36 

std 

0.6291 

0.8667 

2.0131 

1.6039 

1.59% 

10.164 

mean 

0,369 

0.8831 

-0.033 

0.695 

MMMb- 

13.03 

std 

0.8526 

1.2365 

2.5224 

1.9130 

3.21% 

HO 1 


Table 4. Parameter estimates for B 



e, 

e, 

e, 

e. 



true values 

0.21 

1.30 

0.90 

2.00 

noise level 

S/E (dB) 

mean 

0.207 

1.319 

0.881 

1.974 

IWlylli- 

46.97 

std 

0.0506 

0.0488 

0.0575 

0.2477 

0.16% 

9.100 

mean 

-0.058 

2.195 

0.538 

1.789 

jfigfipg 

29.29 

std 

1.7710 

2.8753 

0.3944 

1.1045 


11.076 

mean 

-0.041 

1.159 

-0.020 

1.509 


19.36 

std 

2.0996 

2.6653 

3.6776 

2.0131 


10.164 

mean 

1.349 

-0.076 



Big 

13.03 

std 

3.3518 

4.2351 

38.3601 

2.5224 

3.21% 

5.725 


Table 5. Parameter estimates for 0 


Awilt limulaUd output* end id**l output 


STD of cumulated outputs in different noise levels 




Fig. 16. Model responses 


Fig. 17. Variance in model responses 
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6. Concluding Remarks 

The adaptive weighted least squares algorithm AWLS/MFT has been developed to esti- 
mate the parameters of multivariable linear differential equation models. In addition to es- 
timating the parameters in a least squares setting, it is aimed at finding the best weighting ma- 
trix assuming white measurement noises. As such, it can be regarded as an approximation to 
a maximum likelihood or Gauss-Markov estimate. However, unlike other time-domain based 
approaches, it is frequency-domain based - utilizing only the lowest harmonics in the discrete 
Fourier transform of the data and thereby automatically eliminating high frequency noise if 
the resolving frequency is sufficiently small, i.e., the data interval [O.T] sufficiently long. 
Based upon the comparisons made thus far in simulations and in u tiliz ing real data, the algo- 
rithm appears to compare favorably against some well established time-domain based 
identification techniques like the prediction-error method and an output-error/maximum likeli- 
hood method used for lateral and longitudinal dynamic modeling in aircraft. Applications to 
aircraft modeling with real data has not yet been attempted for the frequency transfer function 
estimation problem, nor for time-varying and nonlinear input/output models, but the exten- 
sions discussed in Section 5 suggest that the potential exists for such applications. 
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